//----------------------------------*-C++-*----------------------------------//
/*!
 * \file   testMesh.cc
 * \author Jeremy Roberts
 * \date   Jul 7, 2011
 * \brief  Tests of various meshes.
 * \note   Copyright (C) 2011 Jeremy Roberts. 
 */
//---------------------------------------------------------------------------//
// $Rev::                                               $:Rev of last commit
// $Author::                                            $:Author of last commit
// $Date::                                              $:Date of last commit
//---------------------------------------------------------------------------//

#include <iostream>
#include <iomanip>
#include <cstdio>
#include <cstdlib>

#include "tran/CartMesh.hh"
#include "tran/Traits.hh"
#include "tran/Definitions.hh"
#include "../utilities/Constants.hh"
#include "../utilities/DBC.hh"
#include "../utilities/Soft_Equivalence.hh"
#include "../utilities/SP.hh"
#include "test/Testing.hh"

using namespace std;
using namespace slabtran;
using namespace util;
using namespace testing;

int test1d_1();
int test1d_2();
int test2d_1();
int test3d_1();

// hard coded test data

// I
const int number_coarse_1 = 4;
const int number_fine_1 = 8;
const double edges_1[] =
{ 0.0, 1.0, 2.0, 5.0, 10.0 };
const double f_edges_1[] =
{ 0.0000, 0.5000, 1.0000, 1.5000, 2.0000, 3.5000, 5.0000, 7.5000, 10.0000 };
const int counts_1[] =
{ 2, 2, 2, 2 };
vector<double> edges1(edges_1, edges_1 + number_coarse_1 + 1);
vector<int> counts1(counts_1, counts_1 + number_coarse_1);

// J
const int number_coarse_2 = 3;
const int number_fine_2 = 6;
const double edges_2[] =
{ 0.0, 2.0, 4.0, 5.0 };
const double f_edges_2[] =
{ 0.0000, 2.0000, 2.5000, 3.0000, 3.5000, 4.0000, 5.0000 };
const int counts_2[] =
{ 1, 4, 1 };
vector<double> edges2(edges_2, edges_2 + number_coarse_2 + 1);
vector<int> counts2(counts_2, counts_2 + number_coarse_2);

// K
const int number_coarse_3 = 2;
const int number_fine_3 = 4;
const double edges_3[] =
{ 0.0, 1.0, 20.0 };
const double f_edges_3[] =
{ 0.0000, 1.0000, 7.333333333333333, 13.6666666666667, 20.0000 };
const int counts_3[] =
{ 1, 3 };
vector<double> edges3(edges_3, edges_3 + number_coarse_3 + 1);
vector<int> counts3(counts_3, counts_3 + number_coarse_3);

int main(int argc, char *argv[])
{
  Require(argc==2);
  int test = atoi(argv[1]);
  cout << " ------ testMesh ------ " << endl;
  cout << " test = " << test;

  int number_of_tests = 5;
  int test_value = 0;

  try
  {
    if (test == 1 || test == number_of_tests) test_value += test1d_1();
    if (test == 2 || test == number_of_tests) test_value += test1d_2() * 100;
    if (test == 3 || test == number_of_tests) test_value += test2d_1() * 10000;
    if (test == 4 || test == number_of_tests) test_value += test3d_1() * 1000000;
  }
  catch (std::exception &err)
  {
    cout << "ERROR: While testing testMesh, "
         << err.what()
         << endl;
  }
  catch ( ... )
  {
    cout << "ERROR: While testing testMesh, "
         << "An UNKNOWN exception was thrown."
         << endl;
  }
  if (test > number_of_tests || test < 1)
    cout << " test " << test << " not available. " << endl;
  if (test_value) cout << " failed with value = " << test_value << endl;
  else
    cout << " passed." << endl;
  return test_value;
}

//! make a 1d mesh using the simple constructor
int test1d_1()
{
  double width = 10.0;
  int count = 8;
  Mesh<_1D> mesh(width, count); //  a derived class object

  const double f_edges[] =
  { 0.0000, 1.2500, 2.5000, 3.7500, 5.0000, 6.2500, 7.5000, 8.7500, 10.0000 };
  vector<double> v = mesh.edges();
  TEST(soft_equiv(f_edges[0], v[0]));
  for (int i = 0; i < number_fine_1; ++i)
  {
    double w = f_edges[i + 1] - f_edges[i];
    TEST(soft_equiv(f_edges[i], v[i]));
    TEST(soft_equiv(w, mesh.width(i)));
    TEST(soft_equiv(1.0 / w, mesh.inv_width(i)));
    TEST(soft_equiv(f_edges[i] + 0.5 * w, mesh.center(i)));
    TEST((i == mesh.convert(i)));
  }

  for (int i = 0; i < mesh.number_cells(); i++)
    TEST((0 == mesh.region_index(i)));

  return 0;
}

//! make a 1d mesh using the less simple constructor
int test1d_2()
{
  Base_Mesh::SP_mesh mesh; // from base smart point
  mesh = new Mesh<_1D> (edges1, counts1);

  vector<double> v = mesh->edges();
  TEST(soft_equiv(f_edges_1[0], v[0]));
  for (int i = 0; i < number_fine_1; ++i)
  {
    double w = f_edges_1[i + 1] - f_edges_1[i];
    TEST(soft_equiv(f_edges_1[i], v[i]));
    TEST(soft_equiv(w, mesh->width(i)));
    TEST(soft_equiv(1.0 / w, mesh->inv_width(i)));
    TEST(soft_equiv(f_edges_1[i] + 0.5 * w, mesh->center(i)));
    TEST((i == mesh->convert(i)));
  }

  const int index[] =
  { 0, 0, 1, 1, 2, 2, 3, 3 };
  for (int i = 0; i < mesh->number_cells(); i++)
    TEST((index[i] == mesh->region_index(i)));

  return 0;
}

int test2d_1()
{
  // make a 2d mesh using the less simple constructor
  Mesh<_2D>::Cell_Edges edges2d(edges1, edges2);
  Mesh<_2D>::Count_Vector counts2d(counts1, counts2);
  Mesh<_2D>* mesh; // from a derived class dumb pointer
  mesh = new Mesh<_2D> (edges2d, counts2d);

  vector<double> vx = mesh->edges(I);
  vector<double> vy = mesh->edges(J);
  TEST(soft_equiv(f_edges_1[0], vx[0]));
  for (int i = 0; i < number_fine_1; ++i)
  {
    double w = f_edges_1[i + 1] - f_edges_1[i];
    TEST(soft_equiv(f_edges_1[i], vx[i]));
    TEST(soft_equiv(w, mesh->width(i, I)));
    TEST(soft_equiv(1.0 / w, mesh->inv_width(i, I)));
    TEST(soft_equiv(f_edges_1[i] + 0.5 * w, mesh->center(i, I)));
  }
  TEST(soft_equiv(f_edges_2[0], vy[0]));
  for (int i = 0; i < number_fine_2; ++i)
  {
    double w = f_edges_2[i + 1] - f_edges_2[i];
    TEST(soft_equiv(f_edges_2[i], vy[i]));
    TEST(soft_equiv(w, mesh->width(i, J)))
    TEST(soft_equiv(1.0 / w, mesh->inv_width(i, J)));
    TEST(soft_equiv(f_edges_2[i] + 0.5 * w, mesh->center(i, J)));
  }

  const int index[] =
  { 0, 0, 1, 1, 2, 2, 3, 3, 4, 4, 5, 5, 6, 6, 7, 7, 4, 4, 5, 5, 6, 6, 7, 7, 4,
    4, 5, 5, 6, 6, 7, 7, 4, 4, 5, 5, 6, 6, 7, 7, 8, 8, 9, 9, 10, 10, 11, 11 };

  for (int i = 0; i < mesh->number_cells_direction(I); i++)
  {
    for (int j = 0; j < mesh->number_cells_direction(J); j++)
    {
      int idx = mesh->convert(i, j);
      TEST((idx == i + number_fine_1 * j));
      TEST((index[idx] == mesh->region_index(idx)));
    }
  }

  delete mesh;
  return 0;
}

int test3d_1()
{
  // make a 3d mesh using the less simple constructor
  Mesh<_3D>::Cell_Edges edges3d(edges1, edges2, edges3);
  Mesh<_3D>::Count_Vector counts3d(counts1, counts2, counts3);
  Mesh<_3D>::SP_mesh mesh;
  mesh = new Mesh<_3D> (edges3d, counts3d);

  vector<double> vx = mesh->edges(0);
  vector<double> vy = mesh->edges(1);
  vector<double> vz = mesh->edges(2);
  TEST(soft_equiv(f_edges_1[0], vx[0]));
  for (int i = 0; i < number_fine_1; ++i)
  {
    double w = f_edges_1[i + 1] - f_edges_1[i];
    TEST(soft_equiv(f_edges_1[i], vx[i]));
    TEST(soft_equiv(w, mesh->width(i, I)));
    TEST(soft_equiv(1.0 / w, mesh->inv_width(i, I)));
    TEST(soft_equiv(f_edges_1[i] + 0.5 * w, mesh->center(i, I)));
  }
  TEST(soft_equiv(f_edges_2[0], vy[0]));
  for (int i = 0; i < number_fine_2; ++i)
  {
    double w = f_edges_2[i + 1] - f_edges_2[i];
    TEST(soft_equiv(f_edges_2[i], vy[i]));
    TEST(soft_equiv(w, mesh->width(i, J)));
    TEST(soft_equiv(1.0 / w, mesh->inv_width(i, J)));
    TEST(soft_equiv(f_edges_2[i] + 0.5 * w, mesh->center(i, J)));
  }
  TEST(soft_equiv(f_edges_3[0], vz[0]));
  for (int i = 0; i < number_fine_3; ++i)
  {
    double w = f_edges_3[i + 1] - f_edges_3[i];
    TEST(soft_equiv(f_edges_3[i], vz[i]));
    TEST(soft_equiv(w, mesh->width(i, K)));
    TEST(soft_equiv(1.0 / w, mesh->inv_width(i, K)));
    TEST(soft_equiv(f_edges_3[i] + 0.5 * w, mesh->center(i, K)));
  }
  // size is 192
  const int index[] =
  { 0, 0, 1, 1, 2, 2, 3, 3, 4, 4, 5, 5, 6, 6, 7, 7, 4, 4, 5, 5, 6, 6, 7, 7, 4,
    4, 5, 5, 6, 6, 7, 7, 4, 4, 5, 5, 6, 6, 7, 7, 8, 8, 9, 9, 10, 10, 11, 11,
    12, 12, 13, 13, 14, 14, 15, 15, 16, 16, 17, 17, 18, 18, 19, 19, 16, 16,
    17, 17, 18, 18, 19, 19, 16, 16, 17, 17, 18, 18, 19, 19, 16, 16, 17, 17,
    18, 18, 19, 19, 20, 20, 21, 21, 22, 22, 23, 23, 12, 12, 13, 13, 14, 14,
    15, 15, 16, 16, 17, 17, 18, 18, 19, 19, 16, 16, 17, 17, 18, 18, 19, 19,
    16, 16, 17, 17, 18, 18, 19, 19, 16, 16, 17, 17, 18, 18, 19, 19, 20, 20,
    21, 21, 22, 22, 23, 23, 12, 12, 13, 13, 14, 14, 15, 15, 16, 16, 17, 17,
    18, 18, 19, 19, 16, 16, 17, 17, 18, 18, 19, 19, 16, 16, 17, 17, 18, 18,
    19, 19, 16, 16, 17, 17, 18, 18, 19, 19, 20, 20, 21, 21, 22, 22, 23, 23 };

  for (int i = 0; i < mesh->number_cells_direction(I); i++)
  {
    for (int j = 0; j < mesh->number_cells_direction(J); j++)
    {
      for (int k = 0; k < mesh->number_cells_direction(K); k++)
      {
        int idx = mesh->convert(i, j, k);
        TEST(idx == i + number_fine_1 * j + number_fine_1 * number_fine_2 * k);
        TEST(index[idx] == mesh->region_index(idx));
      }
    }
  }

  // test insertion of material id.  just use the index above.
  vector<int> matid(index, index + 192);
  mesh->set_cell_id(MATERIAL_ID, matid);
  const vector<int> &matidtest = mesh->get_cell_id(MATERIAL_ID);
  TEST(matidtest.size() == mesh->number_cells());
  TEST(matidtest == matid);


  return 0;
}
//---------------------------------------------------------------------------//
//              end of testMesh.cc
//---------------------------------------------------------------------------//
